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^^ . Abstract 

(N 



A general subtraction scheme, STRIPPER (SecToR Improved Phase sPacE for real Radiation), 

J^ I is derived for the evaluation of next-to-next-to-leading order (NNLO) QCD contributions from 

double-real radiation to processes with at least two particles in the final state at leading order. 

l/^ I The result is a Laurent expansion in the parameter of dimensional regularization, the coefficients 

of which should be evaluated by numerical Monte Carlo integration. The two main ideas are 

(-H I a two-level decomposition of the phase space, the second one factorizing the singular limits of 

Oj' amplitudes, and a suitable parameterization of the kinematics allowing for derivation of subtraction 

' I and integrated subtraction terms from eikonal factors and splitting functions without non-trivial 

flj ■ analytic integration. 

^' 

f\\ , 1. Introduction 

> 

t:;;^- ' Compared to the number of phenomenological applications, where NNLO QCD corrections are 

t'^^- , indispensable, the effort put into the construction of general subtraction schemes for real radiation 

Cn ' at this level of perturbation theory is astounding. The main problems encountered are either, that 

the method is not general and requires tedious adaptation to every specific problem, or that there 

VO ' are many highly non-trivial divergent dimensionally regulated integrals to evaluate. 

^p . This state of affairs is to be contrasted with the comforting situation at the NLO level, where 

general solutions have been available for long. The approach of choice is that of Catani and 
Seymour [l|, later extended to massive states [2| and arbitrary polarizations in real radiation [3|. 
In fact, we would encourage the non-expert reader to consult the original paper [l|, since it gives 
a thorough discussion of all aspects of the problem. The present letter assumes such knowledge. 
^ ' The main features of the Catani-Seymour subtraction scheme are the smooth interpolation of 

the subtraction terms between soft and collinear limits and independence from the phase space 
parameterization achieved by a remapping of phase space points onto the reduced phase space with 
one parton less. There is another scheme at NLO derived by Kunszt, Frixione and Signer (FKS) 
[J|, which is vastly different on the conceptual side. Here, the phase space is first decomposed 
into sectors (originally with the help of the jet function; an independent decomposition has been 
proposed in [5|), and then parameterized with energy and angle variables for easy extraction of 
the subtraction terms. Precisely these ideas will turn out to be crucial for the scheme that we 
shall derive. 

Many approaches have been proposed at NNLO. The most successful ones are Sector De- 
composition [6|, lZI, |8| and Antenna Subtraction j9j| . Sector decomposition is conceptually entirely 
different from the NLO methods cited above. The idea is to derive a Laurent expansion for the 
given amplitude, by first ingeniously parameterizing the complete phase space, mapping it onto 
the unit hypercube, and then dividing it into simplexes, in which singularities are factorized. The 
parameterization is adapted to different singularity structures for different diagram classes. A 
detailed description on the particular example of Higgs boson production can be found in [10|. 
The main drawback is that one has to repeat everything for a new problem, which is relatively 
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easy, if it involves the same number of particles in the final state with the same mass distribution, 
but is expected to be a major effort otherwise. Antenna subtraction, on the other hand, uses com- 
plete matrix elements of simpler processes as building blocks for the subtraction. In this way, the 
integration over the unresolved particle phase space can be made largely with multi-loop methods. 
It is this latter part that involves most work, but the result is general and can be applied to other 
processes without modification. The drawback of this approach is the efficiency loss due to the 
fact that azimuthal correlations characteristic of collinear limits are not taken into account, as the 
simplified matrix elements correspond to unpolarized scattering. There are also other methods, 
which are either specific to a class of problems, such as that of [11], which solves the problem for 
the production of colorless states, or still require the integration of the subtraction terms, as for 
example in 112] and |13j . 

The purpose of this letter is to present a new approach, which should provide a method both 
general and simple to derive. To this end, we need to specify, which problems are essentially 
difficult and which are not. Unlike at NLO, there are in fact two different problems involving 
real radiation at NNLO. One involves the emission of one additional parton (in comparison to 
leading order) out of virtual diagrams. Usually called single-real radiation, since only one parton 
can become unresolved, it can be treated by any of the NLO methods. Of course, the subtraction 
terms require a slight modification. In fact, we believe that the FKS approach will provide the 
result with the least effort. By this argument, we shall ignore single-real radiation and concentrate 
on double-real radiation. The latter problem involves two unresolved partons and only tree-level 
matrix elements. The method presented will not depend on the nature of the initial state. Let us, 
however, consider the more difficult case of hadronic collisions. Due to the factorization theorem, 
the result is a convolution of partonic cross sections with Parton Distribution Functions (PDFs). 
As long as the partonic cross section is an ordinary function, this convolution can be considered 
independently (see Section [^?Tt . Indeed, in actual Monte Carlo generators, one first generates the 
fractions of hadron momenta to be assigned to the partons and then works in the center-of-mass 
frame of the partons, multiplying every event by the PDF weight only at the end. This is the 
point of view that we shall adopt as well. To summarize: we will consider the derivation of the 
Laurent expansion of the double-real radiation contributions for fixed initial parton energies. 

The main concept of our approach is to mix some ideas of the FKS NLO subtraction scheme 
with those of Sector Decomposition. We will decompose the phase space in two stages. At the first 
stage, we will divide the problem into triple- and double-collinear sectors. Then, using an energy- 
angle parameterization of the phase space, we will perform an ordinary sector decomposition 
mimicking the physical singular limits. The crucial novelty is that we will show how to obtain 
general subtraction terms from the last sector decomposition. Here, we will use the knowled ge o f 



NNLO singular behavior of QCD amplitudes as studied in [IJ, [15|, [16| and summarized in [17]. 
While we will not give explicit expressions for the subtraction terms, a task impossible in a letter 
due to the multitude of cases, it is easy to rederive them following the description. 

In the next section, we will derive the scheme on the example of massive particle production. 
The reason for this restriction on the final state is that this letter will be followed by a companion 
publication with process specific information and numerical results for our first application: top 
quark pair production. A subsequent section will, however, present the generalization to arbitrary 
final states. 

2. Massive final states at leading order 

2.1. Phase space 

Let us assume that there are only two massless partons in the final state, the other final 
state particles being massive. In case there is only one massive final state, the presence of soft 
singularities leads to a cross section, which is a distribution in the partonic center-of-mass energy 
squared, s. We will not take this possibility into account, since it has already been extensively 
studied for all processes of phenomenological interest and is a special case involving additional 



complications. The cross section will, therefore, be an ordinary function of s. The considered 
process corresponds to the following kinematical configuration 

P1+P2 ^ ki+k2 + qi-\ ^Qn , (1) 

with 

s^{pi+P2)\ pl^pl^kl^kl^Q, q^^^mj^O, ^ = l,...,n, n>2, (2) 

where pi,p2,ki,k2,qi, ■■■,qn are d-dimensional momentum vectors. The rf-dimensional phase space 
can be written as 



/ 



(3) 
The above definition suggests a factorization of the phase space for this problem into the three- 
particle production phase space of the two massless partons together with an object with invariant 
mass Q^ , and a decay phase space of the composite with momentum Q into the massive particles. 
Such a factorization is motivated by the fact that most divergences are due to the vanishing of 
invariants involving only massless states momenta. The divergences not belonging to this class are 
soft and involve the massive states momenta. In this case, the inverse propagators responsible for 
the singularities vanish proportionally to a linear combination of the energy components of ki and 
^2. This will force us to use these energy components as part of the phase space parameterization. 
At this point the phase space can be written as 



d^n+2 = / -7— 



27r 

dl^-^ki d'^-^k2 d'^-^Q 

(27r)''-i2fcO (27r)''-i2fcO {2ttY-^2QO 

d^'-^q^ 



{27rfS<^''\k,+k2 + Q-pi~P2) 



n (2.).-i2go (2^)'^^'^fa + --- + g»-Q)- (4) 

The integrations over Q can be performed by exploiting the (5-function, which leaves 

d'^'^ki d'^-^k2 



Jd^.,2 = / 



{2TTy-'^2k^ (27r)''-i2fc^ 



^$3 /rf$„(Q), (5) 



where 

Q ^Pi+P2-ki-k2 . (6) 

The same result could have been obtained directly (and trivially) from the original expression 
given in Eq. ([3]), but we wish to keep the interpretation of the integral as described above. Notice 
that the parameterization of d^n{Q) will be of no further concern to us, apart from the fact that 
it is a continuous function of Q. Let us mention, however, that a particularly suitable approach is 
to define d^n{Q) in the center-of-mass frame of Q. In this case, suitable integration parameters 
can be chosen to allow the integration over d^n{Q) to be the first in Eq. ([5|). 

At this point, we can already derive the boundaries of the three-particle phase space. From here 
on, we will work in the center-of-mass system of the colliding particles. Momentum conservation 
from Eq. ([6]) implies 

2{pi + p2)iki + k2) - 2kik2 = s - Q\ (7) 



which can be rewritten as 

2v^(fc° + k°) - 2fc°fc°(l - cos Oi2)^s-Q^ , (8) 

where 612 is the angle between the directions of fci and fc2. Let us introduce the variables 
A, /3, 773,^1,^2 through 



Q' 



s - (Er=i ^if ' 



(9) 



/3 = Ji-SiZ^, (10) 

m = i(l-COS0i2), (11) 

fc°2 - ^/^'^i.2- (12) 

Since Q^ > {J27=i ™i) ^"^ have A, (3, r/3 e [0, 1]. Eq. ^ can now be rewritten as 

a + 6 - /3'%ei6 = A . (13) 

This equation shows that each of the variables a, ^ and 773 can vanish independently, which 
means that the soft and collinear limits are independent for any s. This is due to the fact that 
n > 1, which implied A G [0, 1] (otherwise A — 1). Of course, the independence of the limits can 
be understood intuitively by noticing that they all remove one of the massless partons. Only the 
removal of the two massless partons at the same time (double-soft limit) requires flexibility in the 
choice of Q^ , since then all of the initial state energy is transmitted to the Q system. 

To obtain the upper bounds on the variables, we solve Eq. (|T3)) for one of the energies (the 
expressions are symmetric) 

a.2 = T^^^ • (14) 

1 - p mi2.i 

Independently of the other parton, the maximum is obtained when A — 1, i.e. the composite 
system Q is at threshold. Furthermore, the absolute maximum of the energy occurs, when the 
other parton has vanishing energy 

{^ia)max\i2.l=Q = 1 ■ (15) 

This case does not lead to any divergences, since there is no phase space at threshold for the 
massive Q system. Here, the relative angle, 773, has no relevance. For a finite energy of the 
second parton, the maximum is obtained at 773 = 1, which corresponds to the two partons being 
anti-parallel 



i-e 



2,1 



{il,2)max — ^ _ „2f ' (1^) 

Notice, however, that we will not use this bound explicitly, since our parameterizations will first 
specify the angles and only then the energies, and thus the upper bound will rather be given by 
Eq. P^ with A = 1. Finally, let us note that the range of energy integration is split into ^ > ^ 
and S,2 > Ci at 



1 - y/l ~ I3^rj3 
In fact, we can cover the whole energy integration area by two integration regions defined as 



C I -L - V -L ~ /^ 1/3 



a > 6 : a e [0, 1] , 6 e 



''"""'^^'T^^ 



(18) 



and 

1 - 6 



6 > a : 6 e [0, 1] , 6 G 



0, min ^2 



1 - /32%6 

We will define a function parameterizing the upper bound as follows 



(19) 



?...(0 = min(l,i^2_i_). (20) 

Notice that unless the integration over energies is split as above, the parameterization of the 
phase space cannot be symmetric. Of course, the decisive argument for the split has to do with 
singularities, but we will always keep the symmetry of the expressions, as it mimics the symmetry 
of the final state in the most complicated case of gluons. 

2.2. Decomposition according to collinear singularities 

We will now decompose the phase space taking into account the collinear singularities only. 
The soft singularities will be treated in the next step. 

Collinear singularities occur, when one or more of the following conditions are satisfied: A:i||pi^2, 
fe||Pi,2, fell 1^2- Let us introduce a function Oi{k), such that 9i{k) — 1 if fc||pi, and Oi{k) — 
if fc||p2- In analogy, we introduce 02{k), which satisfies the same requirement after swapping pi 
and p2 ■ Both functions are supposed to fulfill 9i + 02 = 1 and vanish in the respective limits fast 
enough to regulate collinear divergences. The simplest construction satisfying these requirements 
is 9i{k) = 9{k ■ pi). We can now introduce a first partition of the phase space 

1 = (0l(fcl)+^2(fcl))(^l(fc2)+02(fc2)) 

= Ol{ki)0,ik2) + e2{ki)02ik2) + 9,{k,)02{k2) + 92{k,)0i{k2) . (21) 

The first two terms allow for singularities depending only on three of the available momenta. For 
example, the first term will generate singularities when: fci||/i;2|bi, fcilbi, felbi or ki\\k2. This is 
the most complicated case and we will treat it first (see next subsection). We will concentrate on 
the {ki,k2,pi} set, since the other one can be obtained by symmetry. 

The third and fourth terms in Eq. (|2ip will allow for singularities, when fci and fc2 will be 
parallel to different momenta, which means that we can treat the problem as if it were an iterated 
next-to-leading order limit. We must, however, be careful with the situation, when neither of the 
momenta is close to its respective limit, but they both tend to each other. In order to separate 
further this singularity, we will introduce a third function 9^{ki,k2), such that 03{ki,k2) = 1, 
when fci||fc2, and we will assume that ^3(^1,^2) cancels all divergences due to ^1^2 being parallel 
topi^2- 

The final decomposition of the phase space will be given as follows 

'^ ----} — _ 

+02[ki)92(k2) J 

tSSSSIl ~ t[h,li] ] do'^ble-collinear sector 

+{0i{ki)02{k2) + 02{ki)0i{k2))03{ki,k2) } single-collinear sector . (22) 

The first two sectors call for different parameterizations of the phase space adapted to the singular 
kinematical configurations. On the other hand, the third one can be treated along with the first 
sector, since the latter allows for the same single-collinear divergence. Within the first two sectors, 
the parameterizations will be related by symmetry with respect to the exchange of pi O p2. 




Figure 1: Parameterization of the momentum vectors. 



2.3. Triple- collinear sector 

We now turn to the case, where the singularities are due to the three momenta /ci, ^2 and 
Pi. This situation is depicted in Fig. [TJ We will use the recursive d-dimensional definition of the 
momentum vectors 

ri^'^^ = {sine ri('^-^\ cos 0) . (23) 

Using rotational invariance, we define our momenta as follows 



,. ^ ^{f(l,0('^-2)^l) 



Pi 



K = ^i3^Ui,o^''-^\sme,,cose,) , 



^2 



-/3^^2(1,0*'' '^^ sin sin 6*2, COS sin 6'2, COS 6*2) . 



(24) 



Notice that the sign of sin0 is of no relevance for our considerations. We can, for example, assume 
that sin0 > and only restore the sign using a reflection transformation at the end. On the other 
hand, ^1^2 G [OiTt] in complete generality. We will further introduce the notation 



'71,2 = 2(1 -cos 6*1, 2) 



Note that the variable rys defined before is now 

1 



V3 



(1 — cos (f> sin 9 1 sin 62 — cos 9i cos 82) 

(1 — cos(0i — 62) + (1 — cos 0) sin ^1 sin 6*2) 



(25) 



(26) 



The fact that the relative angle 612 vanishes only, when 9i = 62 and = is made explicit in the 
last row. The three-particle phase space can now be written as follows 



(27r)'i-i2fcO (27r)'^-i2fcO 

1 .52-2^/38-8. 



8(27r)5-2£r(l - 2e) 

1 |.i |.i 

dryi (771(1 -?7i))^M dri2 [mi^ - m))'' f dcoscj) [l - cos^ (j))-^'' 
Jo J-i 

JJd^id^2ei~''e2-'\ (27) 

where the energy integration range has been described previously, and d — 4 — 2e. 
As far as the amplitude is concerned, the singular propagator denominators are 

-(Pi-fci)' = s/32^1'71, 
-{Pi-k2f = 5/32^2^2, 

{h+k2f = Sf3%^2m, 

-{Pi-h-k2f = 5/32(^1,71 +e2^2-/3'66%), (28) 

where the signs have been chosen such that the expressions are all positive definite. We have 
omitted the propagators of the massive particles, since they have identical structure, but the 
analogues of the 77 variables can never vanish. The first two structures are already fully factorized. 
Problems arise only in the third and fourth cases. 

Let us focus on 773. We already noted before that the presence of a singularity requires that 
there be 771 = 772 (which is equivalent to Oi =62) and = 0. This is an example of a line singularity, 
since in the three-dimensional space spanned by 771 , 772 and cos 4>, the singularity corresponds to a 
straight line on the cos (j) = 1 plane. On the other hand, a suitable parameterization would only 
exploit two variables, in which case the singular limits would be r/i = 0, 772 = and 771 =7/2, 
independently of the third parameter needed to cover the whole phase space. In order to obtain 
such a parameterization, we will perform a variable change. While the choice is not unique, we 
will use a non-linear variable transformation inspired by that used in [lO| in a similar setting 

C ^ 1 (l-cos(gi-g2))(l+cos0) 

^ 2 1-cos(6li-6l2) + (l~cos(/))sin6lisin6'2 

With this variable, the phase space becomes 

* - 8(2,)^r('i-2.) '""^"- /'"':'«'-':»''- 

d7?id?72 (771(1 - 771)) ''(772(1 - 772)) " 



m -mV ^^ 



//dad^Ci-'^Cz^ 



2e 



(30) 



whereas 773 takes the form 



V. ("^'"^^' . (31) 

7/1 + 7/2 - 27^1772 - 2(1 - 2C) V^l(l - 7/1)772(1 - 772) 



It can be demonstrated that besides the integrable singularities in C, present explicitly in Eq. (|30p 
(which can be remapped to improve convergence), no other singularities are introduced into the 
construction. In other words, C, is not substantial for our discussion. 
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Figure 2; Sector decomposition of the pliase space in the triplc-collinear sector. The variable substitutions, which 
map the integration range onto the unit hypercube are specified. Furthermore, ^2 = imax{^2) and the second 
branch starting with the dashed line is symmetric to the first. 



At this point, we are ready to introduce a further decomposition of the phase space that wiU 
factorize the invariants given in Eq. ([28]) and the ris/lrji — 772 1 factor from the phase space measure 
of Eq. pop into a product of the integration variables and a function, which does not vanish in 
any of the singular limits. The decomposition is presented in form of a tree in Fig. [21 For a given 
ordering of the energies of the massless partons, it contains only five sectors, S(,...,S^. In the 
case of gluons (most intensive from the computational point of view) , there is no need to consider 
the other ordering obtained by symmetry, since the matrix element itself is symmetric. The levels 
in the decomposition tree have a clear physical interpretation: 

I) factorization of the soft singularities; 

II, III) factorization of the collinear singularities; 

IV) factorization of the soft-collinear singularities. 

Each sector specifies an ordering of the relevant variables, thus uniquely defining the limits. 
Collinear singularities require two levels of decompositions, since besides defining which of the 
two partons is allowed to become parallel to pi first, we need to single out the possibility that the 
partons become collinear to each other first and only then to pi . This is achieved at level III. An 
explicit check proves that in all sectors, the necessary factorization is reached. At the same time, 
the inverse propagators of the massive states responsible for soft singularities are treated at level 
I. 



Let us denote the integration measure obtained from Eq. (j30p for a given sector, S, by /i^. 
We consider the contribution of the phase space integral on S to an observable, O, defined by a 
measurement function, customarily called jet function^ Fj, acting on the phase space. We have 



Os^U d(^d7^idindiidi2 fis 0i{ki)ei{k2) / d$„(g) Fj 



M, 



(0) 



(32) 



where A/" is a product of the flux, symmetry factors for the final state, and spin and color average 
factors for the initial state. M}^ is the tree-level amplitude with n + 2 particles. Summation over 
final state polarizations has been omitted in the notation (it is unnecessary, since the formalism 
is correct for polarized amplitudes as well). Moreover, as explained in the Introduction, we have 
ignored the convolution with the PDFs. By construction, we can write the following decomposition 



M5 



M 



(0) 



ri+2 



E 

ai>0 



1 



1 



1 



1 



■Hi 



i>i^ a2 — b2& f^as 



% 



^1 



^2 



'b4e 



Ma^, 



(33) 



where both a,; and bi are integers, bi being defined by fis alone. To make the decomposition unique, 
we require first that Mai,...,ai be regular in the limit of any of the rii^2,£,i,2 variables vanishing. 
Furthermore, if a given A4ai,...,a4, is divided by one of the four variables, i.e. ai > for some i, it 
is not allowed to vanish, when this variable tends to zero. We would now like to obtain a Laurent 
expansion in e for Os • The degree of the singularities given by the at is crucial to determine the 
simplest possible construction. Usual general arguments on the IR structure of QCD amplitudes 
should convince us that Oi < 1, i.e. there are only logarithmic singularities. Thanks to sector 
decomposition, this statement can be verified explicitly for any amplitude, as we will shortly see. 
Let us introduce 



Ms = Vi V2 



■b2e 



fl~b3t(-l- 
?1 ?2 



X fJ-S 



M, 



(0) 
n+2 



With this definition, the observable becomes 



Os=M [ d(:d7j,dij2d^,d^2 e,{k,)6,{k2) ^^^^-T^-ri^ /rf$„(g) Fj m^ 

Jo Vi ' ri2 ^i ^2 J 

and the Laurent expansion is obtained by means of the replacement 



il-be 



b 



E 



(6e)" 



e -^ — ' n\ 

ri=0 



ln"(A) 



where A = 771,2,^1,2, and the "+" -distribution is 



dX 



ln"(A) 



A 



/(A)= f^^-^ifw-fm) 

+ Jo '^ 



(34) 



(35) 



(36) 



(37) 



The /(O) term in Eq. p7p is to be viewed as the subtraction term regularizing the amplitude 
squared, whereas the term proportional to 6{X) in Eq. p6|) is the integrated subtraction term. The 
integrated subtraction terms are integrated over the remaining kinematical variables alongside the 
subtracted amplitude, and no attempt is made at analytic results. After evaluating the integral 
over the (5- function, we can even restore the integration sign, since L dX — 1. Notice finally, 
that the integrated subtraction terms need their own subtraction, which is generated by the same 
expression. The result is obtained by systematically expanding the expression Eq. (j35]) (after the 
substitutions), now a product of polynomials in e (we keep only four terms in Eq. ([36|)), in the 
latter variable down to the finite part. This is correct, since the generated integrals are convergent 
by construction. To some extent, viewing this approach as subtraction terms and integrated 
subtraction terms is inappropriate, but we wish to keep the analogy to the traditional approach. 
Up to now, we have worked with an abstract amplitude, which would, however, have to be 
specified once a definite process would be analyzed. At this point we will make a crucial step in the 



construction, namely the transition to process independent subtraction. Thanks to the decompo- 
sition in terms of physical singularities, each subtraction term generated above and corresponding 
to one or more of the variables vanishing can be obtained in full generality by means of the known 



limiting behavior of QCD amplitudes as summarized in [17[. Indeed, Ms will not vanish in the 
limit, if and only if there is a singularity of the amplitude. The value at this point is obtained at 
the singularity and the factorization into the tree-level amplitude squared with partons removed 
times an eikonal factor or splitting function (divided by the singular invariant) applies. Thus the 
subtraction terms will be given by process specific amplitudes squared (possibly color and/or spin 
correlated) taken in their reduced kinematics, multiplied by the decomposed product of the mea- 
sure fis and eikonal factor or splitting function (divided by the singular invariant). Notice that 
the reduced kinematics is obtained automatically due to our construction of the phase space (the 
concept of momentum mapping used in most subtraction schemes is absent here). Moreover, we 
do not need anymore to think of amplitudes, but only of the eikonal factors and splitting functions. 
Using their form from |17|, it is possible to demonstrate explicitly that Ais is indeed always finite, 
when one or more variables tend to zero. This proves that divergences are only logarithmic. In 
summary, the subtraction terms can be determined once and for all. Due to the multitude of cases, 
the full set is substantially larger than in the next-to-leading order case, but can be derived readily 



with the information provided above and in [17[ by means of simple substitutions (no integration 



is involved, just simple algebra). To be specific, one should proceed in two steps: 

1) Substitute Eqs. (p6| and ([37| into Eq. ([35]), while expanding the "-f" -distributions. The 
result contains 16 different Adg objects with vanishing arguments for 771,2,^1,2 in all possible 
combinations. The term with all variables different from zero is the full amplitude squared 
multiplied by the integration measure. 

2) For a given combination of zero arguments of A4si identify the physical limit and take 
the appropriate factorization formula from [17| . Subsequently, calculate the limit using the 
explicit form of the eikonal factor or splitting function. The reduced tree-level amplitude 
from the factorization formula remains unspecified. 

Finally, let us stress that the subtraction terms are local by construction, a crucial feature to 
assure efficiency. The transverse momentum vectors inducing spin correlations and defining the 
collinear limits can be derived from the momentum parameterization of Eq. p4p , and read in the 
single-coUinear limits 

fc'^i = (0,0(''-4), 0,1,0), 

fc'^2 = (o,o(''-4),2VC(i-C),2C-i,o), 

fc'^3 = (0, Ot-^-^) , yr^, a(l - 2rj) ^,-2a^7j{l- tjK) , (38) 

where fc_Li,2 parameterize the transverse direction of fci,2 with respect to pi, whereas fc_L3 param- 
eterizes the transverse direction of ^2 with respect to fci sd, rji — r]2 — rj (original variables before 
sector decomposition), with a — {r/2 — 'rji)/\V2 ^ J?!]! and we have assumed sint/) > 0. In the triple- 
coUinear limit defined by 771,772 ~^ 0, '7i/^2 =const., the normalization of the transverse vectors is 
not arbitrary. Writing 

k^ = l3^^ip'i + \k1^+0{X^), (39) 

k^ = p'^2P'i + Xkl, + OiX^), 

we can keep k'^^ as above, if A = ^I3^£.i x 2^/771 and 

6 %/P 



kl, = | ^^_^y^^^^ X (0,0('^-^),-2a(l-p)yC(r^,2(l+p)C-(l- Vp)',0) , (40) 
where p = 772/771. 
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Figure 3: Sector decomposition of the phase space in the double-coUinear sector. The notation is as in Fig. [2] 

2.4- Double- collinear sector 

The double-coUinear sector can be treated with exactly the same techniques as the triple- 
colhnear sector. For example, the case where the phase space is restricted by 9i{ki)62{k2){l — 
03(^1,^2)), can be obtained with the previous formulae after the replacement 772 —^ 1 — ^2- Of 
course, there are no singularities depending on the value of 773, which means that one could avoid 
the non-linear variable change. Moreover, the levels II and III from the decomposition tree of 
Fig. [5] can be safely removed. In the end there are only two sectors, as given in Fig. |31 

3. Extension to arbitrary final states 

It may seem at first that the case of arbitrary final states will be substantially more involved 
than the one discussed above. However, the only purpose of the restriction to massive states at 
leading order, was that the system gi,...,g„ would not cause any collinear or soft singularities. 
This may be enforced in the most general case by a decomposition of the phase space similar to 
the one described in Section 12.21 The only difference is that we need more than a constraint on 
the relative angles. 

The double-real radiation at next-to-ncxt-to-lcading order is a correction to a leading order 
process with two massless partons less. The remaining partons need to be well separated and have 
energies bounded from below. For a jet observable, for example, there may be different ways to 
enforce this construction. We will assume, however, that we are given a definite set of initial and 
final states (for a jet observable, there would be a sum over all sub-processes). That is, we need to 
identify a posteriori the leading order processes possible for the given final state. This is achieved 
with the help of a jet algorithm together with the requirement that the final state described by 
jets contain at most two states less than the initial parton configuration. 

In order to be able to use the parameterizations and decompositions of the phase space intro- 
duced previously, we need to extend the division of the phase space described in Section 12.21 Let 
X denote the set of initial state massless partons, and T the set of final state massless partons. 
The total number of final states is again assumed to be at least four. We define two types of 
selector functions acting on the phase space: 9ij^k and 9ij^kh both non- negative, with i,i E F and 
i 7^ j, and k,l G I U J^ and k y^ I, k ^ {hj}, I ^ {hj}- They are required to satisfy the following 
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conditions with m & J- and n G J^ U T 



lim 6l.y_A; : 

fc,.llfc„ 


= 


lim 9jA k 


= 


lim Oij^ki -- 
fe™||fc„ 


= 


lim Oij^ki = 


= 



if {m, n} (t {i, j, k} , 

if rn ^ {i,j} and ?77, is a gluon , 

if {m, n} 7^ {i, k} and {m, n} 7^ {j, /} , 

if rn ^ {i,j} and ?7i is a gluon , 



(41) 



and 



>; 


z2 ^'^o^k + 2I ^^-j'*^' 


pairs 


keiLiJ" leiLiJ" 


^J&r 


fc^{jj} liii.oM 



(42) 



The selector functions thus define a partition of the phase space. Moreover, inclusion of Oij^k 
in the phase space integral allows to obtain a Laurent expansion using the triple-collinear sector 
decomposition described in Section 12.31 assuming that ki and kj take on the role of fci and fc2 , 
whereas kk that of pi. Similarly, Oij^ki allows to obtain a result using the double-collinear sector 
decomposition from Section |2 .41 with ki corresponding to p2- 

While there is ample freedom in defining the above selector functions, it is possible to extend 
the NLO construction from (isj to satisfy the above requirements. Let us define 



d. 



Hjk 



\f2EA 
[\V-sJ 




a 

(1 — COS 6 ij) 


im 




f2Ek\] 


a 

[(1- 



(43) 



[(1 - cose',y)(l - cos6'»fc)(l - cos9jk)Y 



El 



with a, /3 > and 6*^ the angle between ki and fcj. Moreover, let 

hi^k = 

f^ij^k — 



E^ + El 



(44) 



E'k 



El + E] 



El' 



with 7 > and k Cz J^, and hi^k — hij^k = 1 if fc G I. It can readily be verified that the functions 

1 hj2A (45) 



with 



V 





yij.k 


V d,,k ' 




0ij,kl 


1 hi^k hjj 
V dik dji 


>; 

pairs 


>; 

keiuj" 


-dijk ^ dik dji 


i.jeJ" 


k<^{t.j} 


mz,jM 



(46) 



satisfy the conditions in Eq. UT] and Eq. |12] 

We are now able to give the worst case scenario for the number of sectors to consider. At the 
first stage there are {np = |J"| and nj = \I\) 



np 



{nj + nF ~ 2) 



(47) 
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triple-coUinear sectors and 



2( -M( "^ + ^--2 ) , (48) 



double-collinear sectors. Counting both stages, the final number is 

2x5x(^2")^"^ + """'^ + '"'"K"0("'^2"~' ) ■ (49) 

4. Concluding remarks 

The subtraction scheme developed in the previous sections guarantees the possibility to au- 
tomatically obtain double-real radiation contributions to any observable at NNLO. It is general, 
just as much as the NLO subtraction schemes, in the sense that the relevant construction is 
independent of the process. 

Notice that we did not consider it important to show that with this construction observ- 
ables would be finite after combination with the remaining contributions (double- virtual and real- 
virtual). On the one hand, the correctness of the Laurent expansion obtained with the present 
method is proven by the correctness of the intermediate steps. The cancellation of the divergences 
is then guaranteed by the Kinoshita-Lee-Nauenberg and by the factorization theorem. On the 
other hand, aiming at a proof of cancellation leads to a substantial complication of the scheme. 
For us, the simplicity of the numerical implementation was of paramount importance. A similar 
philosophy was already present in the sector decomposition method for phase spaces as devised 
in [tI, |8|. As far as the cancellation of the divergences in dimensional regularization is concerned, 
a comment is, however, in order. In principle, consistency is guaranteed by the use of the conven- 



tional dimensional regularization scheme. This is also implied by our use of formulae from 17 1. 
While our first applications will follow this approach, it is certainly advantageous to use mixed 
schemes, which would allow to compute all tree level matrix elements in four dimensions, just as it 
is often done at NLO. Suitable transition formulae, similar to those of [l9|, will need to be derived 
in the future. 

Finally, let us stress that our scheme is constructed in such a way that the momenta of all 
partons are available together with the integration weight, which allows to obtain arbitrary distri- 
butions on the fly. One may wonder, if the number of sectors will not be a limiting factor to the 
practical feasibility of the calculation. This should not be the case, since one should consider the 
sectors as the analogues of the usual phase space channels. In a Monte Carlo implementation one 
would first choose a sector at random and then the configuration within the sector. The only true 
measure of complexity is the structure of the subtraction terms. Since they have been derived here 
using the physical limits it is expected that their number is minimal, and thus the subtraction 
scheme optimal. Of course, extensions are possible, the most trivial being to add a cut-off on the 
subtraction phase space. Only practice can show, how important such features will turn out to 
be. 
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